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We determine the dynamical critical exponent, «, appearing at the Bose glass to superfluid tran¬ 
sition in two dimensions by performing large scale numerical studies of two microscopically different 
quantum models within the universality class; The hard-core boson model and the quantum rotor 
(soft core) model, both subject to strong on-site disorder. By performing many simulations at dif¬ 
ferent system size, L, and inverse temperature, /3, close to the quantum critical point, the position 
of the critical point and the critical exponents, z, v and rj can be determined independently of any 
prior assumptions of the numerical value of a. This is done by a careful scaling analysis close to the 
critical point with a particular focus on the temperature dependence of the scaling functions. For 
the hard-core boson model we find 2 = 1.88(8), n = 0.99(3) and p = —0.16(8) with a critical field of 
he = 4.79(3), while for the quantum rotor model we find a = 1.99(5), u — 1.00(2) and rj = —0.3(1) 
with a critical hopping parameter of te = 0.0760(5). In both cases do we find a correlation length 
exponent consistent with v = 1, saturating the bound u > 2!d as well as a value of 2 significantly 


larger than previous studies, and for the quantum 
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Most familiar quantum critical points (QCP’s) are 
characterized by Lorentz invariance implying a symmetry 
between correlations in space and time and consequently 
also between the respective correlation lengths ^ ~ Cr • In 
turn, this implies that the dynamical critical exponent, 
defined through ~ C*, is simply 2 = 1 although ^ and 
typically differ by an overall pre-factor. Systems for 
which 2 : yf 1 are comparatively less common and possess 
a distinct anisotropy between space and time since ^ and 
are not only different but now also scale differently. 
One model for which it is generally believed that 2 yf 1 is 
the Bose glass to super fluid (BG-SF) transition describ¬ 
ing interacting bosons subject to disorder, the so called 
dirty-boson problem, modelled by the hamiltonian: 

Hi,h — t ^ ^ J- h.C.^ ^ ^ pr^rA' ~ ^ ^ hr (hr 1) . 

r,e r r 

( 1 ) 

Here e = x, y, and bl,br are the boson creation and 
annihilation operators at site r with hr the correspond¬ 
ing number operator. The parameters of the model are 
the hopping constant t, Hubbard repulsion U, and site- 
dependent chemical potential /ir, inducing the disorder. 

Experimental setups emulating dirty boson physics in¬ 
clude optical lattices [1] adsorbed Helium in random me¬ 
dia [5], Josephson-junction arrays [3], thin-hlm supercon¬ 
ductors [1] as well as quantum magnets such as doped- 
DTN [^. For recent reviews see Ref. HIT]. 

The dynamical critical exponent, 2 , appearing at the 
BG-SF transition has proven exceedingly hard to deter¬ 
mine. Initial theoretical work [5] argued that 2 = d in 
any dimension. This has intriguing implications since 
it implies the absence of an upper critical dimension, 
although a scenario involving a discontinuous onset of 
mean field behavior has beeen proposed [9]. Subsequent 
numerical studies [10IH5) of the system in two dimensions 
were consistent with z = d = 2 and the phase-diagram 


rotor model consistent with z = d. 


is well understood m- More recently, the transition 
in three dimensions has been investigated both numer¬ 
ically [Sin] and experimentally [18] . yielding evidence 
for 2 = d = 3. However, the majority of these numeri¬ 
cal studies were not unbiased since implicit assumptions 
about 2 are needed to fix the aspect ratio L^/(3 in the 
simulations. 

The arguments leading to the equality z = d starts 
with hyperscaling m which states that the singular part 
of the free energy inside a correlation volume is a univer¬ 
sal dimensionless number, {fs/K)Sf‘^T = A. With ^ ^ 
it follows that fs ~ with a finite-size form: [8] 

/.(<5,L,^)~d‘'(‘'+^)E(e/L,^.//3). (2) 

Imposing a a phase gradient dej) along one of the spatial 
directions will then give rise to a free energy difference 
Afs/h = ^p{d(f))^ where p is the stiffness (superfluid den¬ 
sity). Since A/s must obey a form similar to Eq. ([^ 
and since d(j) has dimension of inverse length implying 
dej) ^ 1/^, it follows that p ~ ^ ^u{d- 2 +z) ^ with 

a finite-size scaling form of: 

p = 13/L^), (3) 

Following Ref. [51 an analoguous argument imposing 
a twist in the temporal direction leads to Afs/h = 
^K{dr'p)'^ and k ^ ^ assuming dr4> 

1/^T- On the other hand, it S = (p — pc), the singu¬ 
lar part of the compressibility Ks must from Eq. ([^ obey 
Ks ~ §‘'{d+z)- 2 ^ Assuming that 2 > d, so that k diverges, 
it follows that k = Kg which leads to zv = 1. With z > d 
this result would then contradict the (quantum) Harris 
inequality v > 2/d |2()] invalidating the initial assump¬ 
tion. Hence, one must have z < d. Finally, if one argues 
that for the disordered system n cannot vanish at criti¬ 
cality the relation n ^ S'^ld-z) jjj^plies z = d. 
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FIG. 1. (Color online) Scaling collapse of 142 independent 
simnlations of = jSp for the QR model, (a) Unsealed data 
of versus t for L = 12,..., 32. (b) Scaling collapse of 
the data of panel (a). The data is plotted against the scaling 
variable X = \n{al3/L^)+b{t-tc)L^^''-cL/l3^^^-d{L/l3^^^f 
obtaining a = 1.99(5), u = 1.00(2), tc = 0.0760(5). 


More recent theoretical work dill] has questioned the 
arguments leading to the equality z = d. In particular, 
in the presence of disorder breaking particle-hole symme¬ 
try, it was argued [H] that k — Kg is dominated by the 
analytical background and dr4> ^ 1/Ct should not apply, 
invalidating the relation n ~ leaving z uncon¬ 

strained. A recent numerical [22] study of the hard-core 
version of Eq. ([^ finds z = 1.4 ± 0.05, by analyzing scal¬ 
ing behavior of quantities relatively far from criticality. 
However, in that study, relatively few disorder realiza¬ 
tions (10 — 10^) were used and the location of the QCP 
was not reliably determined. In more recent work, Meier 
et al [23] performed a state of the art calculation of a 
soft-core version of Eq. 0 hnding a signihcantly larger 
value of z = 1.75 ± 0.05 and iz = 1.15(3). While the 
location of the QCP was determined with an impres¬ 
sive precision this latter study does not employ a fully 
quantum mechanical model but instead uses an effective 
classical model for which the temperature dependence, 
at the heart of the scaling with z, is only approximately 
accounted for. It is for instance not possible to calculate 
the specific heat of the underlying quantum model using 
the representation of |23|. 

At present, the value of z at the dirty-boson QCP along 
with many of the other exponents most notably v can 
therefore best be regarded as ill determined, at least for 
the fully quantum mechanical model. It is not known to 
what extent, if any, the relation z = d is violated nor if 
the relation v > 2fd [20| is obeyed, although, as outline 
above it seems reasonable to expect z < d if indeed the 
inequality > 2/d is obeyed and, as shown in |S] ,z>l 
in any dimension. Here we try to answer these ques¬ 



FIG. 2. (Golor online) The Binder cumulant By^r 2 for the QR 
model versus t with /? = Z/^/4. (a) Unsealed data showing a 
crossing close to the critical point G = 0.0760(5). (b) Scaling 
plot versus {t — tc)L^^'' obtained by fitting the data in (a) to 
the form a+b{t—tc)L^^'' +c{t—tc )‘^yielding tc = 0.758(5) 
and V — 0.98(3). 

tions by performing large-scale simulations on two fully 
quantum mechanical models; A hard-core boson model 
(HCB) modelled as a transverse field XY model and a 
soft-core quantum rotor model (QR), both of which are 
subject to strong on-site disorder. In all cases do we use 
10'* — 10® disorder realizations over a large range of tem¬ 
peratures extending down to f3 = 1024 for linear system 
sizes L = 12 — 32. 

Models: The first model we study, closely related to 
Eq. Q, is the quantum rotor (QR) model. It is defined 
in terms of conjugate phase and number operators 9^, 
satisfying [9,.,n,.i] = (5r j.' on a L x L lattice: 

77qr — ^ ^ I COs(0r 0r-|_e) ^ ^ -\ — ^ ^ 71^ (4) 

r,r-|-e r r 

where U is the onsite repulsion, t is the nearest neigh¬ 
bor tunneling amplitude and G represents 

the uniformly distributed on-site disorder in the chemi¬ 
cal potential. As before, e = x, y. The disorder for a 
given disorder realization is not constrained and in all 
simulations we use A = i, 1/ = 1 and tune through the 
BG-SF transtion varying t at constant A. In contrast to 
Eq.0 Tir can take negative as well as positive values and 
one can loosely associtate Ur with deviations from the av¬ 
erage filling rig in Eq. 0> fir — Uq. For convenience we 
study Eq. (|^ using a link-current representation [21] for 
which directed worm algorithms are available [25] . We 
use lattice ranging from L = 12 to L = 32, with 50,000 
disorder realizations for L = 12,..., 28 and 10,000 dis¬ 
order realizations for L = 32 in all cases with 6 x 10^ 
MCS per disorder realization. For the simulations of 
the QR model a temporal discretization of Ar = 0.1 was 
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used, sufRciently small that remaining discretization er¬ 
rors could be neglected. 

The second model we consider is the C/ —>■ oo hard¬ 
core limit of Eq. 0 where the boson occupation number 
is constrained to 0,1. It is therefore therefore equivalent 
to the following S' = 1/2 XF-model on a Lx L lattice in 
a random transverse field: 

H^V = -\Y. + E , (5) 

r.e r 

where € [—h,h] uniformly. In this case we 

tune through the transition by increasing the disorder 
strength, h. Despite the representation as a spin model, 
we shall refer to this as the hard-core boson model 
(HCB). We use a directed loop version of the stochas¬ 
tic series expansion (SSE) [26] to simulate this model. 
This technique does not have discretization errors and 
efficient directed algorithms m HZ] are available. Eor 
the SSE calculations, we use a beta-doubling scheme [55] 
that allows us to very quickly equilibrate at large (3 val¬ 
ues. Eor each temperature in the beta-doubling scheme, 
we average over 48 Monte Carlo sweeps (MCS), with each 
sweep consisting of one diagonal update and Ni directed 
loop updates. Ni is set during the equilibration phase so 
that on average 2{nH) vertices are visited, where nn is 
the number of non-trivial operators in the SSE string. In 
contrast to the QR model, we use here a microcanonical 
ensemble for the disorder by constraining each disorder 
realization to have exactly ^ — 0. This facilitates the 

analysis and is not believed to affect the results |5S| . We 
use at least ^ 10^ disorder realizations per data point, a 
large improvement over |22j . 

In the following [...] denotes the disorder average while 
(...) denotes the thermal average. In simulations of both 
models two independent replicas a, /3 of each disorder re¬ 
alization are simulated in parallel so that combined aver¬ 
ages [(.. .)^] may be correctly estimated as [(.. .)^ (.. .)^]. 

Observables: Our main focus is the scaling behaviour of 
the superfluid stiffness, p, for which the finite-size scaling 
form Eq. 0 was derived. For both models we measure 
p as: 

[{W^ + W^)] 

P= ^ (6) 

where Wx and Wy are the winding numbers in the spatial 
directions. (For the HCB model Eq. 0 is multiplied by 
TT to yield p.) From Eq 0 it follows that Pp = W'^ has 
a particularly attractive scaling form when d = 2, which 
we may write: 

IT" = ^VE((5Li/",L//3i/"), (7) 

where we define 6 = {t — tP) (QR model) and <5 = 
{h — he) (HCB model). We also make extensive use 
of the correlation functions, defined as C(r — r',r — 



O o C{r), L=20, ^ = 800.0 

A A C(r), L=20, ^=800.0 
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FIG. 3. (color online) The correlation functions C{t) and 
C(r) as a function of r, r for a system size L = 20. Results 
are shown for the QR model at the critical point and a range 
of d = 55,..., 800. The solid red line is a fit to /3 = 800 
results for C{t) using the form with 

(z + r])/z=yT = 0.85(2). 


t') = [(exp(f(0r(T) — 6*r'(T'))))] for the QR model and as 
C[t — v',t — t') = [(5'+(r), 5'/7(t'))] for the HCB model. 

Results, QR: A large number of independent sim¬ 
ulations of Eq. Q were carried out at different L = 
12,..., 32 and P = 20,..., 400 close to the QCP. Since 
we expect p to approach zero in an exponential manner 
as L is increased at fixed P and since p is likely expo¬ 
nentially suppressed in the insulating phase it seems rea¬ 
sonable to approximate the function W (x, y) in Eq. 0 
as aexp{f{x,y)) with x = y = LjP^/^. If the 

temperature dependence is carefully mapped out m 
one indeed sees that W{x,y) has a clear exponential de¬ 
pendence. As a first step, we then assume f{x,y) = 
bx — cy — dy^. We can then fit all 142 data points to 
this form determining the coefficients a, b, c, d along with 
te = 0.0760(5) V = 1.00(2) and z = 1.99(5). The results 
are shown in Fig. with a scaling plot using the scal¬ 
ing variable X = ln(a/3/L^) -|- b{t — tP)Lpl'^ — cLjP^!^ — 
d{LlP^I-f. A more refined analysis [301 shows that the 
temperature dependence likely involves a correction term 
= ay^ exp(bx — cy) + dy~'^ exp(6a: — c'y). The correc¬ 
tion term is here proportional to T™ and disappears as 
T tends to zero. It is straight forward to fit all our data 
to this form which yields identical estimates for te-,v, z 
along with w = 0.6(2). Estimating the AIC (Akaike In¬ 
formation Criterion) for the two forms heavily favors the 
latter. 

With a reliable estimate of z we can fix the scaling 
argument L^/P by appropriately selecting /? for each 
L. If we then study the Binder cumulant Hn /2 = 
[(lE‘^)]/[(lT^)]^ we see that at fixed L^ jP it should follow 
a simplified form of Eq. 0: B ^/2 = B{5L^/^). As shown 
in Fig. lines for different L will then cross at tc- This 
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FIG. 4. (Color online) Finite-size scaling analysis of the stiff¬ 
ness for the HCB model using data only at /3 = 512. (a) p 
vs h for system sizes: L — 16, 20,... 32. (b) Crossing of the 
scaling function: R = L^p at the critical point: he = 4.79(3) 
(vertical dashed line). For both panels, the dotted lines are 
second order series expansions of R{x) = a + bx + cx^, fitted 
to the data. 


is indeed the case, confirming our previous estimates. 

Our results for the correlation functions for the QR 
models are shown in Fig. [^for a L = 20 lattice at tc for 
a range of temperatures. Asymptotically, one expects [5] 
C{t) ~ and C(r) ~ r-^d+ 2 -z+rj) _ clearly, 

C(r) drops off much faster than C{t) confirming that 
z 1. However, pronounced finite temperature effects 
are visible in C'(r) arising because the limit /3 ^ has 
not yet been reached and we have not been able to reliable 
determine the power-law describing C (r) . However, from 
C{t) we determine (z -I- ri)/z=yT = 0.85(2) and hence 
rj = —0.3(1) using our previous estimate z = 1.99(5). 

For the QR model we have also verified that the 
compressibility, k, remains finite and independent of L 
throughout the transition, consistent with z < d. Fur¬ 
thermore, a direct evaluation of directly at tc for 
fixed Ifi, expected from Eq. ([^ to scale as ~ , 

yields u = 0.98(4) consistent with our previous results. 

Results, HCB: Due to the hard-core constraint num¬ 
ber fluctuations are dramatically suppressed in the HCB 
model. Combined with the very effective beta-doubling 
scheme we can effectively reach much lower temperatures 
with the HCB model relative to the QR model. Hence, 
we use a simplified form of Eq. ([^ : 

p = ( 8 ) 

suppressing the temperature dependence. We have ex¬ 
tensively verified that this is permissible for the system 
sizes used [30] and that our data appear independent of 
temperature at /3 = 512 to within numerical precision. 
We then ht our data for p at 13 = 512 to an expansion of 


t- 


t- 
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FIG. 5. (color online). The equal-time spatial correlation 
function: C(r) of a 20x 20 lattice for the HCB model at he. (a) 
Convergence of C(r) using the beta-doubling procedure. For 
L — 20, C(r) appears independent of ft for ft > 128. (b) C(r) 
simulation data (open circles) at ft = 512 fitted to the form 
affijr^^ -I- 1/(1/ — r)^'') (dashed line), yielding z + rj =yr = 
1.718(1). 


R in Eq. ([^ to second order: R = a-\-h5L^/'' + c{5L^/’')'^ 
obtaining the estimates: he = 4.79(3), z = 1.88(8), 
v = 0.99(3). The result of this fit is shown in Fig. 
where our simulation data (unfilled markers) is superim¬ 
posed with the quadratic fit (dotted lines). In panel b of 
Fig.H we show the crossing of the scaling function p, at 
he = 4.79(3). By including all results for P < 512 it is 
also possible to perform an identical analysis to the one 
performed for the QR model in Fig.j^ Such an analysis 
similar results for z, tz and he- 

As was the case for the QR model the correlation 
functions show a pronounced temperature dependence 
as shown in Fig. [^a) for C'(r). However, as we lower 
the temperature, C'(r) reaches a stable power-law form 
at ft = 512 for all L studied showing that the regime 
/3 ^ is reached for all L and confirming that the 
temperature dependence can be neglected in Eq. ^ . To 
determine the anomalous dimension p we then fit the re¬ 
sults in Fig. [^a) for L = 20, /? = 512 he = 4.79(3) to 
a power-law form with z -|- p=yr = 1.718(1) as shown 
in Fig. [Kb). Using our earlier estimate of z, we obtain 
p = —0.16(8) in reasonable agreement with the result for 
the QR model. 

For the HCB model we have also calculated the com¬ 
pressibility, K. It remains roughly constant and indepen¬ 
dent of L through the transition. 

Conclusion: Our results for i/ for both models studied 
indicate clearly that v > 2/d is satisfied as an equality. 
For the dynamical critical exponent z, describing the BG- 
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SF transition, we find a value that is significantly larger 
than previous estimates. While there is a slight disagree¬ 
ment in the estimate of z for the two models we studied 
it seems possible that indeed z = d. In light of this it 
now seems particularly interesting to focus attention on 
the transition in d = 4. 

During the final stages of writing this manuscript we 
became aware of Ref. |3T] which for the HCB model reach 
conclusions similar to ours. 
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FIG. 6. Demonstration of the convergence of the spin stiff¬ 
ness: p versus /3 using the beta-doubling procedure, for a 
32 X 32 lattice. As (3 increases, thermal effects can be clearly 
seen to disappear, thereby justifying the use of the simplified 
scaling form in Eq. 



FIG. 7. Gonvergence of spin stiffness: p with number of dis¬ 
order realizations. Nr, for h values used. The shaded region 
indicates the errorbars of points which decrease as more sam¬ 
ples are used. Note that for the equilibration step size chosen 
in our simulations, ~ 10® disorder realizations are necessary 
to obtain true convergence. 


Supplementary material 


Temperature dependence of IF^ = Pp 


HCB simulation details 


A central assumption made in our analysis of the HCB 
results was that a simplified scaling function, Eq. (|^, for 
p. In the main text, we showed the convergence of the 
equal-time spatial correlation function: C(r) in Fig. 
at the critical point using data from the beta-doubling 
procedure. For completeness, we show the convergence 
of the stiffness in Fig. |^for a 32 x 32 lattice. Unlike the 
20 X 20 system, we need to go up to ,0 = 512 before finite 
temperature effects are eliminated. The data indicates 
that for the simulation of even larger system sizes, one 
may need to go up to /3 = 1024 to validate the use of the 
simplihed scaling form. 

Another important point that is often overlooked, is 
the convergence of observable data with the number of 
disorder realizations, Nr used. For the SSE parameters 
chosen, we find that it was necessary to use at least 
~ 5 X 10^ disorder samples before disorder fluctuations 
are reasonably small. This is demonstrated in Fig. [7] for 
the largest lattice size: 32 x 32. For reliable data with 
controlled errorbars, all the HCB data points were av¬ 
eraged over least 10® independent disorder realizations. 
This is contrast to earlier studies [52] most relevant to 
our work, where only 10^ — 10® disorder realizations were 
used. We note that it is quite unlikely that self-averaging 
applies in this model and increasing the lattice size does 
therefore not decrease the number of disorder realizations 
needed. 


Insight into the temperature dependence of can 
be gained by hrst studying p for the QR model without 
disorder. A model for which it is know that z = 1. Re¬ 
sults for p for a 40 X 40 lattice are shown in Fig. [^ for 
/I = 9,... 400. For P L we expect p to go zero in an 
exponential manner while for ^ L it should approach 
a constant. The simplest ansatz is therefore: 


P = 



(9) 


where we have tentatively included an L dependence. 
However, as is clearly evident in Fig. [^ p has a maxi¬ 
mum close to L//3 = 0.4. The presence of this maximum 
signals that there are likely two contributions to p de¬ 
scribing the /3 L and (3 ^ L regimes. (Although we 
note that the existence of two terms does not imply a 
maximum.) We therefore assume the presence of a cor¬ 
rection term proportional to with y > 0. Such a term 
will therefore disappear in the zero temperature limit. 
For our final ansatz we therefore take: 

( 10 ) 

This is the form used in Fig. [^ and it gives an essentially 
perfect fit over the entire range of the figure. We can 
immediately generalize to a scaling form for W'^ = j3p: 

( 11 ) 


= aje 
Jj 


/S.-.4 


mU 


A fit to this form yields w = 0.97(9) in relative good 
agreement with similar correction terms used in Ref. [24j . 
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FIG. 8. (Color online) p as a function of L//3 for the QR 
model without disorder. Results are shown for L = 40. The 
solid yack line indic^es a fit to the L = 40 data of the form 
yielding y — 1.97(9). The dotted line 
indicates the first part of this ht while the dashed line shows 
the second part. 

We now turn to the QR model in the presence of dis¬ 
order. In this case we assume that the scaling variable 
L//3 generalizes to Lj. We therefore expect to find: 

In Fig. results are shown for W'^ for two different lat¬ 
tice sizes L = 12, 20 Assuming z = 2 in this case we plot 
the results against demonstrating that the results 

fall on a single curve with only slight deviations from a 
straight line. It is perhaps surprising that it is the vari¬ 
able that appears as opposed to j{i but this can 

very clearly be verified from the simulations. Perform¬ 
ing a fit to the ansatz we find exceedingly good agree¬ 
ment with the expected form with a correction exponent 
vj = 0.92(7), close to the value for the model without 
disorder. An inspection of our results for p for this sys¬ 
tem size shows that in this case there is no maximum in 
p versus /3. 

We have performed a similar analysis of as a func¬ 
tion of P for the QR model again clearly confirming the 
overall exponential dependence and the presence of the 
two terms. 

Crossing of Pp for the QR model. 

With the z and tf. determined from the ht in Fig. [^/3p 
plotted for different L at a hxed aspect ratio /3 = L^/4 


should also cross at the critical point tc = 0.0760(5). This 
is indeed the case and is shown Fig . [T0| Since this data 
is essentially already shown in Fig IRa) we have in the 
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FIG. 9. (Color online) Pp as a function of Lj^/p for the 
QR model. Results are shown for L = 12,20. The solid 
black line indicates a fit to the L = 12 data of the form 
^ ^ yielding w = 0.92(7). The dotted 

line indicates the first part of this fit while the dashed line 
shows the second part. 



FIG. 10. (Color online) Pp as a function of t for the QR 
model. All simulations are performed using the fixed aspect 
ratio P = L^/4 with z = 2. Lines cross at the critical point 
G = 0.0760(5). 


main text opted to show the crossing using the related 
quantity Byy 2 = [(lF‘*)]/[(bF^)]^ shown in Fig. 

























